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Abstract 

Many endemic species present disjunct geographical distribution; therefore, they 
are suitable models to test hypotheses about the ecological and evolutionary 
mechanisms involved in the origin of disjunct distributions in these habitats. 
We studied the genetic structure and phylogeography of Tibouchina papyrus 
(Melastomataceae), endemic to rocky savannas in Central Brazil, to test hypothesis 
of vicariance and dispersal in the origin of the disjunct geographical distribution. 
We sampled 474 individuals from the three localities where the species is reported: 
Serra dos Pirineus, Serra Dourada, and Serra de Natividade. Analyses were based 
on the polymorphisms at cpDNA and on nuclear microsatellite loci. To test for 
vicariance and dispersal we constructed a median-joining network and performed 
an analysis of molecular variance (AM OVA). We also tested population bottleneck 
and estimated demographic parameters and time to most recent common ancestor 
(TMRCA) using coalescent analyses. A remarkable differentiation among popula- 
tions was found. No significant effect of population expansion was detected and 
coalescent analyses showed a negligible gene flow among populations and an an- 
cient coalescence time for chloroplast genome. Our results support that the disjunct 
distribution of T. papyrus may represent a climatic relict. With an estimated TM- 
RCA dated from ~836.491 ± 107.515 kyr BP (before present), we hypothesized 
that the disjunct distribution may be the outcome of bidirectional expansion of the 
geographical distribution favored by the drier and colder conditions that prevailed 
in much of Brazil during the Pre-IUinoian glaciation, followed by the retraction as 
the climate became warmer and moister. 



Introduction 

The understanding of the ecological and evolutionary mecha- 
nisms involved in the origin of disjunct distributions in plant 
species has been improved by phylogeographic analyses (e.g., 
Azuma et al. 2001; Karanth 2003; Gaudeul 2006; Collevatti 
et al. 2009a). While for some species the disjunct distribu- 
tion represents a climatic relict of the ice ages, caused by 
the range contraction in an ancient more widely distribution 
(e.g., Stehlik et al. 2000; Collevatti et al. 2009a; Lorenz-Lemke 
et al. 2010), long distance dispersal to new suitable habitats 



may also be responsible for disjunct distributions (e.g., Dick 
et al. 2003, 2007; Gaudeul 2006). Also, a number of stud- 
ies based on the fossil pollen record now available for the 
Neotropics show that the distribution and composition of 
vegetation were deeply influenced by the climatic oscillations 
of the Tertiary and Pleistocene. 

During the early Holocene period [until ca. 6000- 
5000 14 C year BP (before present)], the climate was drier in 
most of the South American savannas and distribution of 
savanna-like vegetation in Central and Southeast Brazil was 
more extensive in early compared with the late Holocene 
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(Salgado-Laboriau et al. 1998; Behling and Hooghiemstra 
2000; Behling 2003). In southeastern Brazil, drier climates 
lasted until 8000-8500 14 C year BP or later (Behling 1995, 
2003) and savanna landscapes were dominated by grasslands 
and frequent fires were recorded. The current cerrado vege- 
tation exist in the region only in the latest Holocene period 
(since 970 or 600 year BP for some regions) under the current 
wet climatic conditions, with an annual dry season of about 
4 months. The arid period seems to be a widespread event 
in South America, although the precise age may differ due to 
latitude (Ledru 1993; Salgado-Laboriau et al. 1998; Behling 
2003). Hence, the fossil record shows that savanna expansion 
during glacial periods of the Quaternary was characterized 
mainly by an open savanna with species more adapted to drier 
and highly seasonal climate. In such a dynamic scenario, de- 
mographic expansion and retraction of populations will lead 
to loss in genetic variation and increase in homozygosity, be- 
cause only some genotypes may survive and expand in each 
new cycle (Hewitt 1996; Davis and Shaw 2001). Also, signal 
of bottleneck events may be found in extant populations due 
to retraction during interglacial periods, such as unimodal 
mismatching distribution (Excoffier 2004). Because expecta- 
tion to time to most recent common ancestor (TMRCA) is 
directly related to number of genes in population and prob- 
ability of coalescence is inversely related to the number of 
genes in population (Kingman 1982a, b), most coalescences 
will occur before demographic expansion, when population 
size was small (see Excoffier et al. 2009 for a review). Thus, for 
savanna species adapt to drier and colder climates of glacial 
periods, it is expected that TMRCA coincides with interglacial 
periods when populations became restricted. 

Although a number of studies on phylogeography of 
Neotropical plant species are now available to clarify how cli- 
matic oscillations in Tertiary/Quaternary may have affected 
species distribution in the Brazilian savannas, these studies 
focus mainly on widely distributed species (e.g., Collevatti 
et al. 2003; Ramos et al. 2007). Only one study was per- 
formed with an endemic species restricted to rocky savannas 
from the Brazilian Cerrado Biome (Collevatti et al. 2009a). 
Many endemic species from rocky savannas present disjunct 
geographical distribution and like so are suitable models to 
test hypotheses about the ecological and evolutionary mecha- 
nisms involved in the origin of disjunct distributions in these 
habitats. 

Rocky savannas (cerrado rupestre) are particular vegeta- 
tion communities associated to outcrops of sandstone and 
quartzite soils of the Brazilian Shield in the Cerrado Biome 
in Central and Southeast Brazil, with high number of en- 
demic species and high influence of fire during the dry sea- 
son (Ribeiro and Walter 1998). It is found mainly in Serra 
do Espinhaco, Southeast and Northeast Brazil, and in some 
highlands in Central Brazil such as on Serra dos Pirineus, 
Chapada dos Veadeiros, Serra Dourada, and Serra dos Cristais 



in Goias (GO), in Serra Geral in Tocantins (TO), Chapada 
dos Guimaraes in Mato Grosso, MT, and in similar scattered 
habitats in Chapada da Pratinha (GO and Federal District 
(for detailed description, see Furley and Ratter 1988). 

Tibouchina papyrus (POHL) Toledo (Melastomataceae) 
is an endemic tree species from outcrop quartzite and 
sandstone of the rocky savannas (cerrado rupestre) of the 
Cerrado biome. Flowers are purple and buzz-pollination is 
performed mainly by large bees, such as Xylocopa spp., Bom- 
bus spp., and Centris spp. (Montoro and Santos 2007) and 
the very small seeds are dispersed by autochory. It is con- 
sidered a vulnerable species in the IUCN Red List (IUCN 
2001) and also by the Brazilian Environment Ministry (see 
www.biodiversitas.org.br/ florabr/lista_florabr.pdf ) due to its 
natural endemism and habitat vulnerability. The species is re- 
ported only in three localities in the Brazilian Shield, Central 
Brazil (Fig. 1 ) in Serra dos Pirineus and Serra Dourada, Goias 
(GO), and in Serra da Natividade, Tocantins (TO) (Teixeira 
1969; Montoro and Santos 2007). 

In the present work, we aimed to understand the origin 
of the disjunct geographical distribution of T. papyrus. For 
this, we used a phylogeographical approach to test the hy- 
pothesis that the present disjunct geographical distribution 
is a climatic relict of the ice ages of the Tertiary/Quaternary, 
due to range contraction of a previously more widely dis- 
tributed species caused by the changes in climatic conditions 
that affected suitable habitat distribution. If our hypothesis 
holds, we predicted a high differentiation among popula- 
tions at both chloroplast genome and microsatellite loci, no 
recent gene flow among populations, an ancient TMRCA and 
a signal of bottleneck during interglacial periods. The anal- 
yses were based on the polymorphism at three noncoding 
chloroplast DNA regions and on 10 nuclear microsatellite 
loci. For plants, the comparative analysis of nuclear and or- 
ganelle genomes, with different modes of inheritance and 
mutation and evolutionary rates, may clarify the relative im- 
portance of pollen and seed flow on population structure 
and provide additional insights about the evolution and his- 
torical spread of populations (McCauley 1995; Schaal et al. 
1998; Collevatti et al. 2003, 2009a, b). Additionally, because 
of the haploid nature and mode of inheritance, the effective 
population size of the chloroplast genome is one-half the size 
of the nuclear genome, leading to a stronger effect of genetic 
drift on population genetic structure (Ennos 1994). 

Material and Methods 

Populations, sampling, and DNA extraction 

We sampled T. papyrus populations in each of the three local- 
ities where the species is reported: Serra dos Pirineus (PIR) 
and Serra Dourada (SDO), Goias (GO), and in Serra da 
Natividade (NAT), Tocantins (TO) (Fig. 1). Notice that these 
small number of known local populations is not due to lack of 
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knowledge (i.e., the Wallacean shortfall), because extensive 
sampling of Cerrado plants and rocky savannas were per- 
formed in the last years, searching in many cases especially 
for this species. Because T. papyrus is locally distributed in 
well-delimited patches, we sampled all individuals visualized 
in PIR (216) and SDO (66). In population NAT, we ran- 
domly sampled 192 individuals because this population was 
the largest, with ~1500 individual trees. For each individual, 
expanded leaves were sampled and stored at -80° C. DNA 
extraction followed the standard CTAB procedure (Doyle 
and Doyle 1987). Distances between pairs of populations 
were (Fig. 1): NAT-PIR, 554.82 km; NAT-SDO, 472.32 km; 
PIR-SDO, 146.97 km. 

Sequencing analyses 

For sequencing analyses, we randomly chose 32 adult indi- 
viduals in each population. Three noncoding regions of the 
chloroplast genome (cpDNA) were sequenced: the region be- 
tween the genes psbA and trnH ( Azuma et al. 200 1 ) , t rnS and 
trnG (Shaw et al. 2005), and trnC and ycf6 (Demesure et al. 
1996). Fragments were amplified by polymerase chain reac- 
tion (PCR) in a 20 uL volume containing 1.0 /iM of each 
primer, one unit Taq DNA polymerase (Phoneutria, Belo 
Horizonte, BR), 250 /iM of each dNTP, lx reaction buffer 



(10 mM Tris-HCl, pH 8.3, 50 mM KC1, 1.5 mM MgCl 2 ), 
250 /xg of BSA, and 6.0 ng of template DNA. Amplifica- 
tions were performed using a GeneAmp PCR System 9700 
(Applied Biosystems, Foster City, California, USA) with the 
following conditions: 96°C for 2 min (1 cycle); 94°C for 
1 min, annealing temperature for 1 min (55°C for psbA-trnH, 
56°C for trnS-trnG, 60°C for trnC-ycf6), 72°C for 1 min 
(30 cycles); and 72°C for 10 min (1 cycle). Polymerase chain 
reaction products were sequenced on an ABI 3100 automated 
DNA sequencer (Applied Biosystems, CA) using the DYE- 
namic™ ET terminator cycle sequencing kit (GE Health- 
Care, Uppsala, Sweden), according to manufacturer's instruc- 
tions. All fragments were sequenced in forward and reverse 
directions. 

Sequences were aligned using CLUSTALX (Thompson 
et al. 1997), and characters (each base pair) were equally 
weighted before analysis. For statistical analyses, the data of 
all sequenced regions were concatenated. 

Populations were characterized for nucleotide (jt) and 
haplotype (h) diversity (Nei 1987). Parameters were esti- 
mated using the software Arlequin v. 3.11 (Excoffier et al. 
2005). To understand the origin of disjunct distribution and 
test the hypothesis of vicariance, we performed the following 
statistical analyses. 
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First, we generated a hypothesis of phylogenetic relation- 
ship among haplotypes using the median-joining network 
analysis based on parsimony criteria (Bandelt et al. 1999), 
implemented in the software Network 4.2.0.1 (Forster et al. 
2004). Then, an analysis of molecular variance (AMOVA, 
Excoffier et al. 1992) was performed using Arlequin, to test 
the hypothesis of genetic differentiation among populations. 
The parameter _F S t and pairwise F ST between all pairs of sub- 
population was estimated from AMOVA. Significance of Fst 
and pairwise F S j were tested by a nonparametric permuta- 
tion test (Excoffier et al. 1992) implemented in the Arlequin 
software. Then, the hypothesis that the current pattern of 
haplotype diversity and distribution was caused by contrac- 
tion of an ancient widely distributed population was tested 
under the assumption of a bottleneck followed by a sudden 
expansion. Bottleneck analyses were performed under two 
assumptions: (1) overall populations, under the assumption 
that T. papyrus was more widely distributed in Central Brazil 
during glaciations, and that the present disjunct distribution 
is derived from the contraction of an ancient population (cli- 
matic relict); (2) for each population, based on network re- 
lationships and high differentiation among populations (see 
results below) under the assumption that, after the species has 
attained the disjunction distribution, populations have un- 
dergone independent demographic histories, without gene 
flow among them. The mismatch distribution was obtained 
and the hypothesis of expansion was tested using the Harp- 
ending's Raggeness Index R (Roger and Harpending 1992; 
Schneider and Excoffier 1999) and Tajima's D test of neutral- 
ity (Tajima 1989) also using Arlequin software. 

Finally, a coalescent model (Kingman 1982a, b) was used 
to estimate demographic parameters and TMRCA. The de- 
mographic parameters 9 = 2/iN e (coalescence force or muta- 
tion parameter), g (growth force or exponential growth rate), 
where 9 t = # now exp(-gf/x) and t is time in mutational unit, 
and M = 2N e ml9 (migration force or immigration rate) were 
estimated based on Bayesian estimation using Markov Chain 
Monte Carlo approach implemented in LAMARC 2.1.5 soft- 
ware (Kuhner et al. 2006). Four independent analyses were 
run with 10 initial chains and two final chains and com- 
bined results were generated using Tracer v. 1.5 (Rambaut 
and Drummond 2007). Most probable estimates (MPE) were 
obtained, that is, the highest point on the posterior proba- 
bility curve for a given parameter, which is the best solution 
found by a Bayesian run, and also the confidence interval 
around the estimate of each parameter (Kuhner and Smith 
2007). Time to most recent common ancestor (TMRCA) was 
estimated based on Bayesian phylogenetic analysis imple- 
mented in BEAST 1.6.1 software (Drummond and Rambaut 
2007). A relaxed molecular clock was assumed (uncorrelated 
lognormal), and constant population size (based on the re- 
sults on population expansion, see below). Four independent 
analyses were run for 10 8 generations under HKY+G sub- 



stitution model. Convergence and stationarity were checked 
using the software Tracer vl.5 (Rambaut and Drummond 

2007) . To estimate TMRCA in years, we used mutation rates 
previously estimated for chloroplast noncoding regions (Ya- 
mane et al. 2006). Yamane and coworkers estimated an indel 
mutation rate of 0.8 xlO -9 ± 0.04 x 10 -9 per nucleotide 
per year and 1.52 x 10~ 9 ± 0.06 x 10~ 9 for substitution 
mutation rate. As the chloroplast regions used in the present 
work may have both indels and substitutions we estimated 
TMRCA using each rate and then combined results using 
Tracer v. 1.5. The substitution rates selected here are conser- 
vative, lower than the fastest observed rates for the chloroplast 
genome (e.g., Wolfe et al. 1987). Higher rates would imply 
younger dates for the splitting of the lineages but also the pub- 
lished cpDNA mutation rates are based on grass and herbs, 
which may present twofold faster rates of substitution than 
woody plants (e.g., Clegg et al. 1994; Smith and Donoghue 

2008) leading to an overestimation of divergence times. How- 
ever, we acknowledge a potential limitation of the method, 
which may lead to overestimation of the splitting dates be- 
cause the most recent common ancestor (MRCA) of the hap- 
lotypes (their coalescent) does not necessarily correspond 
to the real temporal split of the populations but may pre- 
cede the actual divergence of the populations (Arbogast et al. 
2002). 

Microsatellite analyses 

For microsatellite analysis, all individuals in each locality 
(474) were genotyped for 10 microsatellite loci, previously 
developed for Tibouchina papyrus (Telles et al. 2011). Micro- 
satellite amplifications were performed in a 15 uL volume 
containing 3.96 /xM of each primer, one unit Taq DNA 
polymerase (Phoneutria, BR), 250 fiM of each dNTP, lx 
reaction buffer (10 mM Tris-HCl, pH 8.3, 50 mM KC1, 
1.5 mM MgCl 2 ), and ~12.5 ng of template DNA. Amplifi- 
cations were performed using a PE 9700 Thermal Controller 
(Applied Biosystems, CA) with the following conditions: 
94°C for 5 min (1 cycle), 94°C for 1 min, 54°C-66°C for 
1 min (according to each locus, see Telles et al. 20 1 1 ) , 72°C for 
1 min (30 cycles); and 72°C for 7 min (1 cycle). The amplified 
products were separated on 6% denaturing polyacrylamide 
gels stained with silver nitrate (Creste et al. 2001) and sized 
by comparison to a 10-bp DNA ladder standard (Invitrogen). 
All individuals were genotyped at least two times in indepen- 
dent PCR amplifications and polyacrylamide gels to avoid 
genotyping error. Although no information on chromosome 
number or ploidy level are currently available for T. papyrus, 
and a high variation has been described for the genus (e.g., 
In = 18 ~ 2n = 56, see Molero et al. 2006), microsatellite 
loci segregation followed the pattern expected for diploid 
species. Thus, the analyses were not affected by the ploidy 
level. 
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Populations were characterized for number of alleles per 
locus, allelic richness based on rarefaction method (see 
Hurlbert 1971; Mousadik and Petit 1996; Petit et al. 
1998), and observed and expected heterozygosities un- 
der Hardy-Weinberg equilibrium (Nei 1978), and inbreed- 
ing coefficients if). Analyses were performed with FS- 
TAT 2.9.3.2 (Goudet 2002) and randomization based tests 
with Bonferroni correction were performed for deviation 
from Hardy-Weinberg expectations and linkage equilibrium 
(Goudet etal. 1996). 

Genetic differentiation among populations was assessed 
by Wright's F-statistics, -Fit, -Fst. and F\ S (Wright 1951), 
obtained from an analysis of variance of allele frequencies 
(Cockerham 1969; Weir and Cockerham 1984), implemented 
in FSTAT 2.9.3.2 (Goudet 2002). To verify the contribution of 
stepwise-like mutations to the genetic differentiation, we esti- 
mated Slatkin's R S j (Slatkin 1995) obtained from an analysis 
of variance of allele size following Goodman (1997). Then, 
we tested the hypothesis that F S t = Rsj based on Hardy et al. 
(2003), using the software SPAGeDI (Hardy and Vekemans 
2002). The comparison of _F S t and R S j provides insights into 
the role of drift and mutation in population differentiation, 
because _R S t is expected to be larger than F S t under stepwise- 
like mutations, but equal when differentiation is caused only 
by drift (see Hardy et al. 2003). 

The effect of strong reduction in population size due to 
bottleneck in current population diversity was analyzed us- 
ing the Wilcoxon sign-rank test (Luikart et al. 1998) imple- 
mented in Bottleneck software (Cornuet and Luikart 1997). 
The bottleneck process may cause a faster loss of heterozy- 
gosity under mutation-drift equilibrium (H eq ) than loss of 
observed heterozygosity (H 0 ). Hence, populations that have 
experienced recent reduction in effective population size may 
present higher H 0 than H eq for a given number of alleles in 
the population (Maruyama and Fuerst 1985). The analyses 
were run under the assumption of the three models of mu- 
tation to verify the sensitivity of the results to the mutation 
model (see Cornuet and Luikart 1997): IAM (Infinite Allele 
Model), SMM (Stepwise Mutation Model), and TPM (Two- 
phase Mutation Model). 



Demographic parameters were estimated for microsatellite 
data based on the coalescent model as described above for se- 
quencing data, correcting for diploid genome (0 = 4fiN e and 
M = AN c ml6), using LAMARC software. Time to most recent 
common ancestor (TMRCA) was estimated from overall 9. 
We used the lowest mutation rate reported for microsatellite 
markers in plants, 2.4 x 10 -4 mutation per allele per gener- 
ation (95% CI = [1.4 x 10" 4 ; 4.2 x 10~ 4 ]) (Thuillet et al. 
2002), often quoted in the range of 10 -3 to 10 -4 per locus 
per generation (e.g., Ellegren 2000; Udupa and Baum 2001; 
Thuillet et al. 2002; Vigouroux et al. 2002; Marriage et al. 
2009). We chose the lowest value based on the comparison of 
-Fst and Rst (see results below) and used a generation time 
of eight years (MPC Telles, unpublished data). 

Results 

Sequencing analyses 

Amplification of the noncoding regions psbA-trnH, trnS- 
trnG, and trnC-ycf6 generated fragments of 268 bp, 
580 bp, and 294 bp, respectively. Al sequences were submit- 
ted to GenBank database (accession numbers: JN969604 to 
JN969699, for psbA-trnH, JN969700 to JN969795, for trnC- 
ycf6; JN969796 to JN969892, for trnS-trnG). 

For the 96-sequenced individuals, psbA-trnH noncoding 
region presented eight different haplotypes, with five poly- 
morphic sites. For the region trnC-ycf6, we found three hap- 
lotypes and two polymorphic loci. For trnS-trnG, we found 
five haplotypes and four polymorphic sites (Table SI). Gene 
diversity for psbA-trnH (h = 0.771) was very similar to the 
value found for trnS-trnG (h = 0.775) and nucleotide diver- 
sity {it = 0.0046; SD = 0.0033) was higher than for trnS-trnG 
(jt = 0.0025; SD = 0.0017). The region trnC-ycf6 presented 
the lowest level of gene diversity (h = 0.341) and nucleotide 
diversity {it = 0.0023; SD = 0.0019). The combined data from 
the three sequence regions generated a fragment of 1 142 bp, 
with 17 haplotypes, and 11 polymorphic sites (Table SI). 

Population SDO presented the highest genetic diversity 
(Table 1), followed by PIR. However, nucleotide diversity 



Table 1. Genetic diversity and demographic parameters for Tibouchina papyrus populations based on Bayesian coalescence analysis for combined 
chloroplast DNA data from psbA-trnH, trnS-trnG, and trnC-ycf6. N = 32 for all populations (Pop), fc-number of haplotypes; ft-haplotype diversity; 
jr-nucleotide diversity; SD-standard deviation; 9-coalescent parameter; SE-standard error; g-exponential growth parameter (values followed by ns 
did not differ from zero based on the confidence interval around the estimate). 



Pop 


Combined cpDNA 












k 


h 


jt (±SD) 




9 (SB) 


g(SE) 




NAT 


4 


0.591 


0.0013 ± 


0.0009 


7.83E-05 (7.39E-06) 


309.839" 


' (46.7671) 


PIR 


6 


0.762 


0.0012 ± 


0.0009 


5.09E-05 (3.24E-06) 


322.301" 


' (100.1479) 


SDO 


7 


0.837 


0.0012 ± 


0.0008 


1.14E-04 (7.54E-06) 


353.308" 


'(113.5003) 


Overall 


17 


0.912 


0.0029 ± 


0.0017 


2.77E-04(1.16E-06) 


394.452" 


' (57.9844) 
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Figure 2. (A) Median-joining network based on the sequences of chloroplast noncoding regions from 96 individuals of Tibouchina papyrus from 
three localities in Central Brazil. Circumference size is proportional to the haplotype frequency. Haplotypes corresponds to those described in S1 . All 
mutations are shown in the network, mv1 , median vector. Different colors were assigned for each population: NAT, black; PIR, red; SDO, blue. (B) The 
model of bidirectional expansion of ancestral population from PIR to SDO and NAT. The arrows show the possible paths following the hills with two 
possible paths from PIR to NAT. 



was very similar among populations (Table 1). Populations 
did not share any haplotype (Fig. 2). Haplotype H07 from 
population PIR was the most ancient (Fig. 2), but H02 was 
the most frequent haplotype. 

Analysis of molecular variance showed a high differen- 
tiation among populations (F ST = 0.684, P < 0.001), and 
pairwise F$j were significant and high for all pairs of popu- 
lations (NAT-PIR, F ST = 0.821, P < 0.001; NAT-SDO, F S t 
= 0.802, P < 0.001; PIR-SDO, F ST = 0.755, P < 0.001). 

Roger-Harpending mismatch distribution test and 
Tajima's neutrality test did not show any significant effect 
of population expansion after a bottleneck (P > 0.05 for 
all populations and overall populations). Coalescent analysis 
performed with LAMARC software also showed that popu- 
lations are not expanding but have constant size with very 
low mutation parameter 9 (Table 1), indicating a historical 
low effective population size or a low mutation rate. Num- 
ber of migrants per generation was negligible for all pairwise 
comparisons (Table 2). 

Coalescent analyses performed with BEAST software indi- 
cated an ancient TMRCA for T. papyrus, dated from 836.491 
year BP ± 107.515 year BP for chloroplast genome (Fig. 3). 



Table 2. Number of migrants per generation for Tibouchina papyrus 
populations, based on Bayesian coalescence analysis for combined 
cpDNA data from psbA-trnH, trnS-trnG, and trnC-ycf6, and on 1 0 micro- 
satellite loci. Migration direction is from populations (Pop) in the lines 
into populations in the columns. 



Pop 


Combined cpDNA 




Microsate 


ites 




NAT 


PIR 


SDO 


NAT 


PIR 


SDO 


NAT 




0.0137 


0.0091 




0.0048 


0.0007 


PIR 


0.0184 




0.0079 


0.0574 




0.0041 


SDO 


0.0076 


0.0054 




0.0014 


0.0007 





The ancestral diverged in two major lineages that originated 
population SDO and the ancestral of PIR and NAT popu- 
lations (Fig. 3). The last lineage, diverged at 661,342.29 ± 
100.515 year BP. 

Microsatellite analyses 

All pairs of microsatellite loci were in linkage equilibrium (all 
P > 0.05). Genetic diversity was higher at population NAT, 
but allelic richness was higher in SDO (Table 3). Populations 
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Figure 3. Relationships and TMRCA (time to most recent common ancestor) of the lineages from 96 individuals of three populations of Tibouchina 
papyrus in Central Brazil, based on Bayesian coalescent analyses of chloroplast sequences. Haplotypes corresponds to those described in Table SI. 
Different colors were assigned for each population: NAT, black; PIR, red; SDO, blue. Time scale is in years before present. 



Table 3. Genetic diversity and demographic parameters for Tibouchina papyrus, based on the polymorphism at 10 microsatellite loci. W-number 
of individuals analyzed; A r — allelic richness; H e -expected heterozygosity; H D -observed heterozygosity; /-inbreeding coefficient (values followed by ns 
are not significant, for P= 0.00167, Bonferroni correction for a nominal value of 0.05); 6-coalescent parameter; SE-standard error; g-exponential 
growth parameter (values followed by ns did not differ from zero based on the confidence interval around the estimate). 



Population 


N 


A, 


He 


H 0 


f 


0(SE) 


g(SE) 


NAT 


192 


1.8 


0.409 


0.386 


0.057™ 


0.015 (0.0003) 


999.896(1.1061) 


PIR 


216 


1.5 


0.205 


0.155 


0.246 


0.009 (0.0005) 


999.981 (0.1034) 


SDO 


66 


2.1 


0.357 


0.280 


0.216 


0.019 (0.0009) 


-56.648 ns (5.4947) 


Overall 


474 










1.505 (0.0007) 


5.433 (1.3657) 



PIR and SDO presented observed heterozygosities statisti- 
cally lower than expected under Hardy-Weinberg equilib- 
rium (Table 3). 

Populations were highly differentiated (P ST = 0.712, P = 
0.002). R ST also showed a high level of genetic differentia- 
tion (Rst = 0.831, P < 0.001), but f S T and R S t were not 



significantly different (P = 0.427). A significant amount of 
inbreeding was found (Pjs = 0.127, P = 0.002), and non- 
random mating among populations (Fu = 0.748, P = 0.002). 
Pairwise Fst were also significant and high for all pairs of 
populations (NAT-PIR, F ST = 0.843, P < 0.017; NAT-SDO, 
Fst = 0.691, P < 0.017; PIR-SDO, F ST = 0.655, P < 0.017). 
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Bottleneck analysis showed an excess of H 0 in relation to 
H eq , which is a signal of population expansion after a drastic 
size reduction for the three mutation model (P < 0.001). 
Coalescent analysis also showed population expansion but 
not for population SDO (Table 3). Mutation parameter (9) 
was low in all populations and migration among popula- 
tions was negligible (Table 2). Using the coalescent parame- 
ter overall population (0 = 1.505), we were able to estimate 
the TMRCA, which dated to 50,166.67 year BP (95% CI = 
28,666.67-86,000.00 year BP). 

Discussion 

Our results showed that the disjunction distribution of 
T. papyrus is most likely due to the range contraction of 
an ancient more widely distribution, representing a climatic 
relict of drier and colder ice ages in Central and Southeast 
Brazil. 

Populations of T. papyrus are highly differentiated for both, 
chloroplast and microsatellite markers. High differentiation 
is usually found in endemic species (e.g., Collevatti et al. 
2009a; Zhou et al. 2010; Ge et al. 2011; Jones and Gibson 
2011; Wang et al. 2011). Pollination and dispersal systems of 
T. papyrus may also foster population isolation and differ- 
entiation. This species is pollinated by large bees, and seeds 
are very small and dispersed by autochory. Although large 
bees may fly long distances (Goulson and Stout 2001), the 
patchy distribution of T. papyrus may favor the isolation of 
pollinators on patches of plants, restricting pollen flow due to 
long distance among suitable habitats. Also, the small seeds 
of T. papyrus are not winged and might not fly long distances, 
which may also pose limitations on the expansion of popu- 
lations by diffuse or jump dispersal. Thus, pollen dispersal 
most likely contributes more to long distance dispersal than 
seeds. Dispersal may be constrained because the species is un- 
able to cross barriers or because species is habitat specialist 
and may not succeed in habitats available in the barrier (Cox 
and Moore 2005). Tibouchina papyrus life history suggests 
that both factors may constrain species expansion. Hence, a 
modification in habitat availability should be necessary for 
species geographical distribution expansion. 

Despite the low mutation parameter 0, populations pre- 
sented high haplotype diversity for chloroplast genome. This 
is most likely due to the historical constant population size 
with no signal of severe reduction in population size, which 
maintained diversity despite low mutation rate. However, 
analyses of microsatellite loci showed a different demograph- 
ical history and lower diversity than that usually found for 
microsatellite loci in savanna trees (e.g., Collevattti et al. 2001, 
Braga et al. 2007; Pereira et al. 2008; Rabelo et al. 2011). In 
fact, analyses on demographic history for microsatellite loci 
showed a reduction in population size that was not detected 
by chloroplast genome. This is most likely due to the differ- 



ences in mutation rate and the time since the bottleneck event 
occurred (see Excoffier et al. 2009 for a review). Because mi- 
crosatellite loci display high mutation rates, typically in the 
range of 10 -3 to 10 -4 per locus per generation (see Goldstein 
and Schlotterer 1999 for a review), they are more prone to the 
effects of recent demographic events. Those differences are 
denoted in the TMRCA, which was lower for microsatellite 
loci than for chloroplast genome. Besides differences in mu- 
tation rate, this may also be the outcome of a more recent bot- 
tleneck that may not be detected by a region with slower evo- 
lutionary rate. As coalescence of microsatellites which evolves 
faster coincides with ~50 kyr BP (~28-86 kyr BP), that pre- 
dates the last glacial maximum (LGM), we hypothesize that 
populations suffered a bottleneck during the Sangamonian 
interglacial period (~1 10-130 kyr BP) and started to expand 
during the last glaciation (~ 12-1 10 kyr BP). Also, because the 
power to detect expansion after a reduction in population size 
is dependent on time since population reduction when us- 
ing neutrality test or mismatch distribution (Excoffier 2004; 
Excoffier et al. 2009), ancient demographic expansion may 
be undetectable in both chloroplast genome and micro- 
satellite loci. However, major coalescence times (see Fig. 3) 
coincided with pre-IUinoian period, specially the Cromerian 
interglacial (~455-620 kyr BP) which may indicate periods 
of smaller effective population size. 

Network analysis indicated that the haplotype H7 (Fig. 2) 
from Serra dos Pirineus population is the extant more an- 
cient haplotype and despite the distance between Serra dos 
Pirineus and Serra Dourada, the ancestor of the last popula- 
tion diverged first (Fig. 3). This result may be the outcome 
of limited genome sampling leading to errors in coalescence 
time estimation and should be interpreted with conscious. 
However, despite the difference between the median-joining 
network based on parsimony criterion and on Bayesian coa- 
lescent tree, the clade comprising the lineages of PIR popula- 
tion has the lower number of mutations and thus the lowest 
number of character changing, suggesting that PIR is the an- 
cestral population of T. papyrus. Hence, we hypothesize a 
bidirectional population expansion from Serra dos Pirineus 
toward the Southeast and Northwest during the glacial pe- 
riods. This region of Central Brazil is comprised by many 
highlands, such as Chapada Pratinha, Chapada dos Vead- 
eiros that could have connected PIR population to Serra da 
Natividade, the Northern distribution of T. papyrus. 

Although estimations presented high variation, the TM- 
RCA for populations dates from the Pre-IUinoian Cromerian 
interglacial that occurred between 620,000 and 455,000 year 
BP. This interglacial was followed by the Anglian Glaciation 
(Pre-Illinoian), that last from 455,000 to 300,000/380,000 
year BP (Gibbard and Kolfschoten 2004; Gibbard et al. 2007). 
In the last glacial period, the Wisconsin Glaciation (~1 10,000 
to 12,000 year BP), the expansion and changing in geo- 
graphical distribution is observed in the fossil records for 
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some species of high altitudes, such as Podocarpus (e.g., 
Colinvaux et al. 1996), Aracauria, and other conifer taxa 
(Ledru 1993). Based on the above scenario, we hypothesize 
that T. papyrus expanded its distribution in the glacial peri- 
ods of the Pre-Illinoian or Cromerian complex, favored by 
the colder and drier conditions that created suitable habitats. 
With the warmer and moisture conditions of the following 
interglacials, many populations extinguished and T. papyrus 
became restricted to drier quartizite and sandstone soils, at- 
taining the present distribution in the Central Brazil. Thus, 
the present disjunct distribution may represent a climatic 
relict of an ancient widely distributed population. Indeed, 
the high differentiation and low gene flow between popula- 
tions are also evidences that vicariance is responsible for the 
disjunct distribution. This same process was also suggested 
for Lychnophora ericoides (Asteraceae), a species from rocky 
savannas that co-occurs with T. papyrus in SDO and PIR 
(Collevatti et al. 2009a), despite the more recent TMRCA of 
L. ericoides. The differences in TMRCA are most likely due 
to the more wide geographical distribution of L. ericoides 
in the rocky savannas of the Central Brazil and differences 
in demographic parameters and ongoing gene flow among 
populations (Collevatti et al. 2009a). 

In conclusion, our results strongly support that the disjunct 
distribution of T. papyrus may represent a climatic relict and 
that long distance gene flow is unlikely. Tibouchina papyrus 
has undergone a recent bottleneck that may have caused loss 
in genetic diversity. Since migration among populations is 
negligible, the probability of local extinction due to demo- 
graphic and genetic stocasticity may be very high. Tibouchina 
papyrus is a rare and endemic species of sandstone outcrop 
habitat, which is highly unstable and suffers high levels of dis- 
turbance caused by fire during the dry season and sandstone 
and quartzite disruption, mainly during the rainy season. 
These disturbances may be highly variable among local habi- 
tats and may cause sudden modifications in population size 
(Collevatti et al. 2010). An important outcome of our study 
is the evidence that T. papyrus populations are three differ- 
ent ESUs (Evolutionary Significant Units, Moritz 1994) that 
may have adaptive differences. Hence, these populations re- 
quire separate genetic management. These results are highly 
important for planning and executing scientifically sound 
conservation programs for the species. Moreover, due to re- 
markable differentiation and ancient divergence of the three 
lineages from the three localities, we suggest a more in depth 
taxonomical study of this species to review its classification in 
only one species. Finally, the present work showed remarkable 
concordance in the mechanisms involved in the evolution of 
the disjunct distribution of two unrelated species from rocky 
savanna, despite differences in some demographic parame- 
ters. This result is a strong evidence for the hypothesis raised 
in this work for the evolution of disjunct distribution in plant 
species from rocky savannas in Central Brazil. 
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